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ABSTRACT 

We construct a secular theory of a coplanar system of A/^-planets not involved in strong mean 
motion resonances, and which are far from collision zones. Besides the point-to-point New- 
tonian mutual interactions, we consider the general relativity corrections to the gravitational 
potential of the star and the innermost planet, and also a modification of this potential by 
the quadrupole moment and tidal distortion of the star. We focus on hierarchical planetary 
systems. After averaging the model Hamiltonian with a simple algorithm making use of very 
basic properties of the Keplerian motion, we obtain analytical secular theory of the high order 
in the semi-major axes ratio. A great precision of the analytic approximation is demonstrated 
with the numerical integrations of the equations of motion. A survey regarding model parame- 
ters (the masses, semi-major axes, spin rate of the star) reveals a rich and non-trivial dynamics 
of the secular system. Our study is focused on its equilibria. Such solutions predicted by the 
classic secular theory, which correspond to aligned (mode I) or anti-aligned (mode II) apsides, 
may be strongly affected by the gravitational corrections. The so called true secular resonance, 
which is a new feature of the classic two-planet problem discovered by Michtchenko & Mal- 
hotra (2004), may appear in other, different regions of the phase space of the generalized 
model. We found bifurcations of mode II which emerge new, yet unknown in the literature, 
secularly unstable equilibria and a complex structure of the phase space. These equilibria 
may imply secularly unstable orbital configurations even for initially moderate eccentricities. 
The point mass gravity corrections can affect the long term- stability in the secular time scale, 
which may directly depend on the age of the host star through its spin rate. We also analyze 
the secular dynamics of the V Andromede system in the realm of the generalized model. Also 
in this case of the three-planet system, new secular equilibria may appear. 

Key words: celestial mechanics - secular dynamics - relativistic effects - quadrupole mo- 
ment - analytical methods - stationary solutions - extrasolar planetary systems 
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1 INTRODUCTION 

The currently known sample of extrasolar planet|^ comprises of 
many multiple-planet configurations. Their architectures are di- 
verse and, usually, very different from the Solar system config- 
uration. Some of them consist of so called hot-Jupiters or hot- 
Neptunes, with semi-major axes ~ 0.05 au and the orbital periods 
of a few days. The multiple-planet systems may also contain com- 
panions in relatively distant orbits. Such systems are non-resonant 
(or hierarchical). 

The discovery of multi-planet systems and their unexpected 
orbital diversity initiated interest in their secular evolution. Gen- 
erally, the analysis of secular interactions, including tidal effects, 
general relativity (GR), and quadrupole moment (QM) corrections 
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to the Newtonian gravitation (NG) of point masses, follows the the- 
ory of compact multiple stellar systems ( M ardling & Lin|2002l|Na-| 
[gasawa & Lin||2003| . The secular evolution of planetary systems 
has been already intensively studied in this context by many au- 
thors. For instance, [Adams & Laughlin] p006b| ) consider the GR 
corrections in the sample of detected extrasolar systems. They have 
shown that the GR interactions can be particularly important for the 
long-term dynamics of short-period planets. Also |Adams & Laugh-| 
[iml ( |2006a| l; |Mardling| ( [20b 7 ) study the evolution of such objects 
taking into account the tidal circularization of their orbits, and sec- 
ular excitation of the eccentricity of the innermost planet by more 
distant companions. These secular effects imply constraints on the 
orbital elements of the innermost orbits, including yet undetected 
Earth-like planets, which are expected to be found in such systems. 
The quadrupole moment and relativistic effects may detectably af- 
fect the light-curves of stars hosting transiting planets in a few years 
time-scale ( .Miralda-Escude 2002, Agol et al.. 20051 |Winn et al.] 
|2005| ). In general, the interplay of apparently subtle perturbations 
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with the point mass Newtonian gravity depends on many physical 
and orbital parameters, and it may lead to non-trivial and interest- 
ing dynamical phenomena. Still, the problem is of particular im- 
portance for studies of the long-term stability of the Solar system 
(see very recent works by |Benitez & Gallardo|2008[[Laskar 2008 ). 

Because a fully general study of such effects would be a very 
difficult task, we introduce some simplification to the planetary 
model. We skip tidal effects caused by dissipative interactions of 
extended star and planetary bodies, in particular the tidal friction 
(TF) damping the eccentricity and the tidal distortion (TD) that 
modifies planetary figures, and could influence their gravitational 
interaction with the parent star. The tidal effects are typically much 
smaller than the leading Newtonian, relativistic and quadrupole 
moment contributions (Mardling {1001) . Also their time- scale is 
usually much longer than of the most significant conservative ef- 
fects. In this work, we focus on the secular dynamics of planetary 
systems over an intermediate time scale relevant to the conserva- 
tive perturbations. Our goal is to obtain a qualitative picture of the 
secular system that may be useful as the first order approximation 
to more general model of the long-term planetary dynamics. 

Moreover, we show in this paper that even if we skip the non- 
conservative tidal effects then the GR and/or QM corrections to 
the point mass Newtonian gravity can be much more significant 
for the secular dynamics than the mutual NO interactions alone. 
These effects may change qualitatively the view of the dynamics 
of planetary systems as predicted by the Laplace-Lagrange theory 
(Murray & Dermott 2000 ) or its recent versions ( e.g.,|Lee & Peale| 
2003 , Michtchenko & Malhotra 2004 , Michtchenko et aL||2006] 
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The model without dissipative effects may be investigated 
with the help of conservative Hamiltonian theory. Assuming that 
planetary orbits are well separated and the system is far from mean 
motion resonances and collision zones, we can apply the averaging 
proposition ( Arnold et al.|1993| ) to derive the long term evolution of 
their orbital elements. This approach can be classified among secu- 
lar theories having a long history in the context of the Solar system 
( [Brouwer & Clemence|196T| [Murray & DermottpOOO] ). 

The present work is a step towards a generalization of the sec- 
ular model of coplanar 2-planet system by Mic htchenko & MaT] 
[hotra[ ( |2b04) and its analytical version app lied to A/^-planet sys- 
tem in ( Migaszewski & Gozdziewski|2008a| ). To explore the phase 
space globally without restrictions on eccentricities, we simplify 
the equations of motion through the averaging of perturbations to 
the Keplerian motion, with the help of semi-numerical method pro- 
posed by Michtchenko & Malhotra, ( ,2004^ ). To obtain analytical re- 
sults, we also use a very simple averaging algorithm that can be 
applied to perturbations dependent on the mutual distance of inter- 
acting bodies (Migaszewski & Gozdziewski 2008a). These works 
demonstrate that the secular evolution can be precisely described 
in wide ranges of the orbital parameters, including eccentricities 
up to 0.8-0.9, for well separated (hierarchical) orbits with small ra- 
tio of the semi-major axes, a ~ 0.1. Here, the secular NG theory is 
very helpful to derive the more general model including the GR and 
QM effects. Basically, without the averaging, the only possibility 
of investigating the long-term secular dynamics relies on numeri- 
cal solutions of the equations of motion ( [Mardling & Lin[[2002[ ). 
However, due to extremely different time scales which are related 
as days (the orbital periods of inner planets) to 10^-10^ years of 
the secular orbital evolution, the numerical integrations are of very 
limited use when we want to investigate large volumes of initial 
conditions rather than isolated orbits. In such a case, the analyti- 



cal theory can be very helpful to get a deep insight into the secular 
dynamics. Moreover, when necessary, particular solutions can be 
studied in detail with the help of the direct numerical integrations. 

The plan of this paper is the following. In Sect. 2 we formu- 
late the generalized model of a coplanar, A/^-planet system, account- 
ing for the relativistic and quadrupole moment corrections to the 
NG-perturbed motion of the innermost planet. To make the paper 
self-consistent, we average out the perturbations to the Keplerian 
motion with the help of the averaging algorithm described in ( |MP] 
gaszewski & Gozdziewski 2008a,) . In Sect. 3 we apply the ana- 
lytical and numerical tools to study the influence of the GR and 
QM corrections on the secular evolution. In particular, we recall 
the concept of the so called representative plane of initial condi- 
tions ( Michtchenko & Mal hotra 2004 ). We focus on the search for 
stationary solutions in the averaged and reduced systems (periodic 
orbits in the full systems), and we investigate their stability and bi- 
furcations with the help of phase diagrams. In this section we also 
derive interesting conclusions on the stability of the unaveraged 
systems. In Sect. 4, we study the phase space of two-planet sys- 
tems in wide ranges of parameters governing their orbital configu- 
rations (masses, semi-major axes ratios, eccentricities) and physi- 
cal parameters (flattening of the star). To illustrate the application 
of the secular theory to multi-planet configurations, we consider 
the three-planet \) Andr system, and we investigate the QM (stellar 
rotation) influence on its secular orbital evolution. 



2 GENERALIZED MODEL OF TV-PLANET SYSTEM 

The dynamics of the planetary system can be modeled by the 
Hamiltonian function written with respect to canonical Poincare 
variables (see, e.g., [Poincare| 1 8971 [Laskar & Robutel|1995| >, and 

expressed by a sum of two terms, 9{ = + i^ert, where 



epl " 



2p; 



(1) 



Stands for integrable part comprising of the direct sum of the rel- 
ative, Keplerian motions of N planets and the host star. Here, the 
dominant point mass of the star is mo, and mi <^ mq, i = I, . . . ,N 
are the point masses of the A/^-planets. For each planet-star pair we 
define the mass parameter jUi = {mQ-\-mi) where k is the Gauss 
gravitational constant, and P/ = (1/m/ + l/mo)~^ are the so called 
reduced masses. We consider the perturbing Hamiltonian i^ert as a 
sum of three terms. 



■^ert — + -^R + 



(2) 



where is for the mutual point-mass interactions between plan- 
ets, is for the general (post-Newtonian) relativity corrections 
to the Newtonian gravity, and takes into account the dynami- 
cal flattening and tidal distortion of the parent star. 

We consider the secular effects of i^o, which can be ex- 
pressed as follows: 



i=l j>i 



iiiiij 



Vi Vj 



A/j- mo 
direct part indirect part 



(3) 



where r/ are for the position vectors of planets relative to the star, 
p/ are for their conjugate momenta relative to the barycenter of 



the whole {N + l)-body system, A/^ ^ 
distance between planets / and j. 



\vi — Vj\\ denote the relative 
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For hierarchical planetary systems (weakly interacting bina- 
ries), the secular time scale of Newtonian point-mass interactions 
may be as short as 10^ years, up to Myrs. In general, these interac- 
tions cause slow circulation of the apsidal lines. The GR correction 
to the NG potential of the star-planet system also leads to the cir- 
culation of pericenters. Moreover, the time scale of this effect may 
be comparable to that one forced by mutual Newtonian interactions 
between planets. In such a situation, one should necessarily include 
the relativistic corrections to the model of motion. 



2.1 General relativity corrections 

The GR correction will be applied only to the innermost planet, 
so we skip the direct GR perturbations on the motion of the more 
distant companions, as well as the mutual GR interactions caused 
by planetary masses. Usually, the GR corrections to the Newtonian 



potential are expressed in terms of the PPN formalism ( e.g.,|Kid- 
jder^ 1 995 .) . A lternatively, we found a very clear paper of |Richar~ 
|son & Kelly] ( |1988j ) conveniently providing explicit Hamiltonian of 
the two body problem with the GR term. Following these authors, 
= ^^r/P (i.e., rescaled by the reduced mass) may be 
written as follows: 



^R = yiP'+y2- 



-73 



(r-pr 



1 



(4) 



where yi , 72 , Ys , 74 are coefficients defined through 



Yi 



(1-3V) 



8c2 ' 2c2 

and c is the velocity of light in a vacuum, /i = k^{mQ +mi), V = 
momi/(mo + mi)^, r is the astrocentric distance, and P is the as- 
trocentric momentum of the innermost planet (normalized through 
the reduced mass): 



A^(3+v) 



73^ 



2c2' 



14- 



1 



279 274 
47i(v-v)v+ -^v+ ^(r-v)r 



(5) 



where v = r stands for the astrocentric velocity of the innermost 
planet (the relativistic corrections from other planetary bodies in 
the system are skipped). Hence, in the relativistic Hamiltonian, we 
put P = V with the accuracy of 0(c~^) and then the Hamiltonian is 
conserved up to the order of 0{c~^). 

2.2 Rotational and tidal distortions of the star 

Fast rotating stars are significantly flattened and that in turn may 
lead to important deviations of the NG potential. In the absence of 
a close planet and non-radial pulsations, the star has the rotational 
symmetry, regarding its mass density and shape. The gravitational 
potential of such a body can be expanded in harmonic series ex- 
pressed in terms of the Stokes coefficients. The well known prop- 
erty of this expansion is that the gravitational potential of rotation- 
ally symmetric bodies retains only terms related to the so called 
zonal harmonics: 

^pm = ^f://f^yP/(sin(^), (6) 

where Rq stands for the characteristic radius of a sphere encom- 
passing the body, and it can be fixed as the equatorial radius of the 
star, P/(sin(|)) is the Legendre polynomial of the l-th order, (|) is the 
astrocentric latitude and // are non-dimensional Stokes coefficients, 
/ > 1 . It can be shown that for rotationally symmetric objects, Ji = 
for / odd. Basically, // (and the leading term with J2, in particular) 



can be determined numerically with the help of the theory of stellar 
interiors combined with a model of rotation and helio- seismic data 
( |Godier & Rozelot|1999[|Pijpers|1998] . To calculate these coeffi- 
cients, one must know the mass density as well as the hydrostatic 
figure of the star. There are attempts to estimate J2 for stars hosting 
planets. For inst ance, [ Iorio|(|2006|) determ ined quadrupole moment 
of HD 209458 ( Charbonneau et al.|2000| ) as J2 ^ 3.5 x 10"^ how- 
ever with large error of ~ 10"-^, making the result not very credible. 
In general, the estimates of // are very uncertain even for the well 
explored Sun. For this slowly rotating dwarf (with rotational period 
of ~ 30 days), J2 ~ 10~^ with an upper bound of 10~^. Moreover, 
in the literature we found approximations of current J2 of the Sun 
in the range up to 10~^. 

Due to indefiniteness of higher order zonal harmonics, we 
consider the dynamical effects of the leading J2 term only. For a 
coplanar system, in which planets move around the star in its equa- 
torial plane, the first non-zero term in harmonic series, Eq. [6] has 
the form of: 

1 ..^91 
'2 



(7) 



The above formulae may be conveniently expressed ( [Mardling &| 
|Lin|2002'] ) in terms of the stellar spin frequency: 



-hRl{l+mi/mo)kLa^\ 



(8) 



where ki is the tidal Love number ^Murray & Dermott|2000^ which 
can be related to J2 through: 



h 

1 



mo 



(9) 



q is the ratio of centrifugal acceleration and gravitation at the sur- 
face of the star, and Q. is the spin rate. Hence, we have: 



J2 



3 k^niQ 



(10) 



In this work we adopt the standard value of = 02 ( [Nagasawa &| 
|Lin|2005| . Obviously, the dynamical flattening is more significant 
for fast rotating stars. Before the Zero Age Main Sequence (ZAMS) 
stage, the rotational periods may be as low as a few days, down to 
1 day for young (~ 100 Myr) Sun-like stars. The rotational periods 
of ~ 40 days are typical for 8 Gyrs old, evolved objects. According 
to the scaling rule of J2^^^, the zonal harmonics of young objects 
may be as large as 10~^. Because the rotational period may change 
by two orders of magnitude during the life-time of the star, also 
the quadrupole moment may change by a few orders of magnitude. 
Hence, flattening caused by fast rotation not only can affect signifi- 
cantly the orbital acceleration which can compete with the GR and 
mutual NG contributions but it may also introduce a dependence of 
the system dynamics on the age of the parent star (Nagasawa & Lin| 
|2005| see also Sect. 5 and Sect. 6). 

We also consider a contribution of the tidal bulge (TB) caused 
by the presence of point-mass innermost planet, assuming that only 
the extended star is distorted due to the tidal interactions. Then the 
Hamiltonian is corrected by the following term ( [Mardling & Lin| 
|2002l ): 



-^Rlk^mi (1 +mi/mo)^L 



1 



(11) 



Both gravitational corrections are then = -^^pin + -^^b- 

We account for direct dynamical effects related to only 
in the motion of the innermost planet although these perturbations 
can be also quite easily included for the remaining star-planet pairs. 
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In the realm of the secular theory of non-resonant, hierarchical sys- 
tems, we assume that other companions are much more distant from 
the star. For the semi-major axes ratio ~ 0. 1, the GR-i-QM perturba- 
tions acting on such outer companions are by orders of magnitude 
smaller than for the innermost planet. However, the dynamics of 
the innermost body influences indirectly the secular dynamics of 
the whole system. This will be demonstrated in this work in the 
cases of two- and three-planet configurations. We also underline 
that with the simplified model of interactions, we can study the dy- 
namics of quite compact star-planets configurations because they 
can be parameterized, in general, by individual planetary masses, 
semi-major axes, and physical parameters of the star. 



2.3 The secular model of non-resonant planetary system 

To apply the canonical perturbation theory, we first transform H to 
the following form: 



i?/(I,(^)=i^epl(I) + ^ert(I,(^), 



(12) 



where (I,(|)) stand for the action-angle variables, and i^ert(I,(|>) ~ 
854epi(I). where 8 <C 1 is a small parameter. In this work, we ap- 
ply the well known approach to analyze the equations of motion 
induced by Hamiltonian, Eq.[T2] that relies on the averaging propo- 
sition (see, e.g., [Arnold et al. 1993 ). By averaging the perturbations 
with respect to the fast angles {the mean longitudes or the mean 
anomalies) over their periods, we obtain the secular Hamiltonian 
which does not depend on these fast angles. Simultaneously, the 
conjugate momenta to the fast angles become integrals of the secu- 
lar problem. In the planetary system with a dominant stellar mass, 
the orbits (apsidal lines) also slowly rotate due to mutual interac- 
tions, hence the longitudes of periastron and the longitudes of node 
become slow angles. Assuming that no strong mean motion res- 
onances are present, and the system is far enough from collisions, 
the averaging makes it possible to reduce the number of the degrees 
of freedom, and to obtain qualitative information on the long-term 
changes of the slowly varying orbital elements (i.e., on the slow 
angles and their conjugate momenta). 

The transformation of the Hamiltonian to the required form 
(Eq.[T2| may be accomplished by expressing this Hamiltonian with 
respect to the (modified) Delaunay canonical elements. These vari- 
ables can be related to the Keplerian canonical elements (Murray] 
|& Dermott,2000| ). Actually, we use the following set of canonical 
action-angle variables which can be obtained after an appropriate 
canonical transformation of the Delaunay elements: 




(13) 



Hi = Gi (cos /, - 1 



where Mi are the mean anomalies, a/ stand for canonical semi- 
major axes, ei are the eccentricities, // denote inclinations, 05/ are 
the longitudes of pericenter, and 12/ denote the longitudes of as- 
cending node. The choice of 05/ instead of CO/ is important for fur- 
ther applications, because ^2/ are undefined (and irrelevant) for the 
dynamics of the coplanar system. We note that the geometrical, 
canonical elements (a/, ^/, //, 05/, 11/) may be derived through the 
formal transformation between classic (astro-centric) Keplerian el- 
ements and the relative Cartesian coordinates (e.g., F erraz-Mello| 
|et al.|2006{|Morbidelli|2002| , with appropriate rescaling of the as- 
trocentric velocities. In the settings adopted here, the Cartesian co- 
ordinates are understood as Poincare coordinates, i.e., astrocentric 



positions of planets, and canonical momenta taken relative to the 
barycenter of the system. 

The A/^-planet Hamiltonian expressed in terms of the modified 
Delaunay variables (T3] ) has the form of: 



• -^jY ^ -^^^ ■>h->Gi, gi , Hi , hi ) 



i=l 



i=l,...,N 



In this Hamiltonian, // play the role of the fast angles. In the absence 
of strong mean motion resonances, these angles can be eliminated 
by the following averaging formulae: 

\ rlTi r2% 

^^^ = 7^^ •••/ ^GdMi...dMN+ (14) 
(27i)^^ Jo Jo 



i=l,...,N 



1 r27i 
'2n 



r2% I r2% 

/ ^^dMi + — / ^udMi. 
Jo 271 Jo 



Hence, we can rewrite the secular Hamiltonian in the symbolic 
form of: 



(15) 



Now, we try to average out each component of this Hamiltonian 
expressed with respect to the fast angles (the mean anomalies). 

2.4 Averaging Newtonian point-to-point interactions 

A simple averaging of the Hamiltonian of the classic planetary 
A'^+ 1-body model is described in our previous paper ( |Migaszewski| 
|& Gozdzi ewski|2008a ). The algorithm makes use on the very ba- 
sic properties of the Keplerian motion and it relies on appropriate 
change of integration variables in Eq. [14] It may be also applied to 
perturbations expressed through powers of the relative distance. 

To recall the main result, the secular Hamiltonian of the A^- 
planet system can be described as a sum of two-body Hamiltonians 
evaluated over all pairs of planets: 



(^g)=II«^^). (16) 

The direct part of the disturbing Hamiltonian reads as follows: 
k^mim j 



/=2 Vl-^} 



^ ^^^'^^(^/,^,-,A05/j 



(17) 



The explicit formulae for functions ^*^''^\^/,^j, A05/j) are given 
in ( Migaszewski & Gozdziewski|2008a| ). It is well known that the 
indirect part averages out to a constant and it does not contribute to 
the secular dynamics ( [Brouwer & Clemence|1961| ). 

2.5 Averaging the PPN relativistic potential 

Making use of the formulae for the relativistic PPN Hamiltonian 
by Richards on & Kell'y] ( |1988| ), we write down the mean relativistic 
potential as follows: 



,4\ /{r-xY 



+ Y4 



(18) 



with the accuracy of 0{c ^). We should average out each com- 
ponent of this Hamihonian over the mean anomalies. It appears 
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to be possible with the algorithm in fMigas zewski & Gozdziewski| 
|2008a| l. To average out the whole PPN Hamiltonian, we must cal- 
culate a few integrals written in the general form of: 



\ r2% I fill 



271 Jo 271 . 

where J is a scaling function defined with: 



1- 



df (1 + ^cos/)^ 



(19) 



(20) 



and / is the true anomaly of the innermost planet (note that we skip 
index "1" of that planet). Components of the mean Hamiltonian can 
be written explicitly as follows: 



r) ^ 

2 



(r-v) 



1 



Functions !fi and are given through: 



r2n 

Jo (T 



sinV 



1 

2nJo (i+ecos/) 
'■2t sin^/ 



df, 



-1 



df. 



271 JO 1+^cos/ 
These integrals can be calculated in the closed form: 



3(2-^2_2y/l_^2 



2^4 



Finally, the secular relativistic potential is the following: 



(21) 
(22) 
(23) 
(24) 

(25) 
(26) 

(27) 
(28) 



Because all secular Hamiltonian corrections we account for in this 
work (see below) do not depend on the mean anomalies, Li 2 are 
constants of motion. It also means that the second term in Eq. |28] 
reduces to the constant and does not contribute to the long-term 
dynamics. Writing down the secular Hamiltonian with respect to 
the action-angle variables, we have: 



- const. 



(29) 



Here, the action variables (L, G) are defined with usual formulae 
known in the Keplerian motion theory. We recall that in the sec- 
ular PPN Hamiltonian, only the star-innermost planet interactions 
are considered. The perturbation has been also scaled by the re- 
duced mass. To properly add the GR star-inner planet interactions 
to the Hamiltonian of A/^-planets, we should account for the mass 
factor. Because the canonical Delaunay elements of the A/^-planet 
system are defined through Eq.[T3] hence, to obtain the proper scal- 
ing, we should multiply the two-body Hamiltonian, Eq.|29] through 
p= l/(l/mo + l/mi). Finally, we obtain the following secular 
Hamiltonian related to the GR correction: 



- const. 



(30) 



In this Hamiltonian all slow angles are cyclic, hence L and G would 
be integrals of motion in the absence of mutual planetary interac- 
tions. However, the element G is no longer constant when we con- 
sider the NG contributions. The evolution of canonical angle 05, 
reads as follows: 



G5gr ■- 



^(^R) _ 3a/4|35 



3A/3/2 



c2l3g2 c2a5/2 (1 



(31) 



This is the well known formulae describing the relativistic advance 
ofpericenter. We derive it here to keep the paper self-consistent. In 
the above formulae, we skip the symbol (...) of the mean, which 
basically should encompass the canonical angle 05 and the conju- 
gate actions (L, G) as well as the symbols of orbital elements (a, e). 
However, we should keep in mind that after calculating the average 
over the fast angles, these actions/elements have the sense of the 
mean actions/elements. 



2.6 The effect of quadrupole moment of the star 

To average out the quadrupole moment of the star, we have to cal- 
culate integrals from the astrocentric distance of the planet taken 
in the power of "-3" for the spin distortion, and the power of "-6" 
for the tidal distortion. To calculate these averages, we express the 
astrocentric distance r of the inner planet with respect to the true 
anomaly and we replace the integration variables dM = Jdf. For 
/ > 1, we find 



, L-i I , V-(cos-/), (32) 



where, for s even, the average of (cos^ /) over the mean anomaly 
reads as follows: 



{cosV> 



2^71 



r(i-|)^r(5+i) 



(33) 



while for s odd, the average over the mean anomaly (cos^ /) = 0. 
The average of the leading term in (i^pin) has the form of: 



PkL{l^mi/mo)RlQ? 

6^3(1-^2)3/2 



(34) 



Having this secular Hamiltonian, we can calculate the apsidal fre- 
quency forced by the spin-induced QM of the star: 



G5st 



dG ^ 2G^L^ 
Similarly, we calculate the leading term in {^jb}' 

^kiRl (1^4 ^ 3^2 ^ 1) k^mi (1 + mi /mo) 



(35) 



x9/2 



(36) 



and the apsidal frequency due to this correction: 

a(i^TB) _ 15pl3y^^/^5 (^G4-14L2G2 + 21L4)miA/^ 



05tb = 



dG 



SGi^L^mo 



(37) 



The tidal bulge induced by a close companion may be very im- 
portant in some configurations. In fact, the TB effect may be even 
dominant over the spin-induced quadrupole moment perturbations. 
For instance, for t) Andr b it is ~ 20 times as big (see |Mardling &| 
Lin2002, Nagasawa & Lin|2005| for details). 
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3 TOOLS TO ANALYZE THE SECULAR DYNAMICS 



3.1 Representative plane of initial conditions and equilibria 



To overview the dynamical influence of the GR and QM corrections 
on the point mass NG dynamics, we compare the secular evolution 
of the eccentricity of the innermost planet in the sample of detected 
multiple extrasolar systems. In this experiment, we basically follow 
[Adams & Laughli n (2006b ), however, solving the averaged equa- 
tions of motion by means of the numerical integrator. Moreover, 
we carry out this test to illustrate the different and rich behaviors 
of multi-planet systems in the detected sample; a systematic anal- 
ysis is described in subsequent sections. The results are shown in 
Fig.[T] Each panel in this figure is labeled with the parent star name. 
The orbital elements of the selected systems, and the parameters of 
their parent stars come from the Jean Schneider Encyclopedia of 
Extrasolar Planet^ In this experiment, we fixed rotational period 
of the star, T,^, = 30 days, ki = 0.02 (hence J2 - lO-^-lO^^). In 
some cases, the differences between predictions of the generalized 
theory and by the classical model are significant, in some other 
cases both theories give practically the same outcome. Looking at 
the elements of the examined systems, we may conclude that the 
corrections become very important for systems with the innermost 
planet close to the star, and with other bodies relatively distant. 

In fact, it is already well known (e.g., [Adams & Laughlin] 
[2006b[ l that the additional effects are important when the strength 
and characteristic time scale of perturbations stemming from the 
GR and QM become comparable with the secular time- scale of mu- 
tual NG interactions. Then the GR and QM regarded as perturbing 
effects may induce significant contribution to the secular evolution 
and the interplay of these effects with the NG interactions may lead 
to very complex and rich dynamics. As we will see below, this con- 
dition is particularly well satisfied for strictly hierarchical systems, 
e.g., HP 21 7107, HD 38529 (Fischer et al.'2001"), HD 190360 fWogi^ 
et al.|2005|), HD 11964 (Butler et al. 2006), v Andromedae (Butler 



etal. 



1999 ). The secular time scales can be also comparable when 



the planetary masses are relatively small because the NG-induced 
apsidal frequency decreases proportionally to the products of these 
masses while, in the first approximation, the GR and spin-induced 
apsidal motion of the innermost orbit does not depend directly on 
the planetary mass (see Eqs. [37][35] ). 

Still, the results illustrated in Fig.[T]have limited significance 
for the study of the global dynamics. Drawing one-dimensional 
time-orbital element plots, we can analyze the dynamical evolu- 
tion only for a few isolated initial conditions. This can be a serious 
drawback, if we recall that the initial conditions of the discovered 
systems are still known with large uncertainties. Some critically im- 
portant parameters governing the secular evolution, like the masses, 
nodal lines and inclinations of orbits are poorly constrained by the 
observations or undetermined at all. There is also a technical prob- 
lem: the direct integrations over the secular time-scale are CPU in- 
tensive. Yet looking at isolated configurations, we obtain only a 
local view of the dynamics. Instead, following the methodology of 
Poincare, we should try to understand globally the perturbing ef- 
fects and the resulting dynamics. We should investigate the whole 
families of solutions rather than a few isolated phase- space trajec- 
tories. In such a case, the application of secular analytical theories 
become critically important. 



The simplest class of solutions which can be studied most effec- 
tively, are the equilibria (or stationary solutions). In the multi- 
parameter dynamical systems, their positions (coordinates), stabil- 
ity, and bifurcations provide information on the general structure of 
the phase space. Hence, we focus on the stationary solutions emerg- 
ing in the secular model of a two-planet system with the GR and 
QM corrections. 

After the averaging of ^ over the mean anomalies (i.e., over 
the orbital periods), the secular Hamiltonian (i^ec) does not de- 
pend on li = Mi anymore. Therefore, as we mentioned already, the 
conjugate actions L/ become constants of motion (hence, the mean 
semi-major axes are also constant). Moreover, in the coplanar prob- 
lem the longitudes of nodes are undefined and irrelevant for the dy- 
namics. The secular energy of two-planet system does not depend 
on individual longitudes of pericenters 05/, but only on the their dif- 
ference AG5 = G5i — 032 ( [Brouwer & Clemence|196l| ). These basic 
facts can be expressed with the help of appropriate canonical trans- 
formation ( Michtchen ko & Malhotra.2004j : 



AG5 = G5i —052, 
^2, 



C = Gi+G2. 



(38) 
(39) 



Because 052 is the cyclic angle, the conjugate momentum equal to 
the total angular momentum C is conserved. For a fixed angular mo- 
mentum as a constant parameter, the phase-space of reduced system 
become two-dimensional, and the system is integrable. We can now 
choose A05 as the canonical angle and Gi as the conjugate momen- 
tum. Alternatively, the role of the momentum may be attributed to 
ei . Obviously, G2 (or ^2) becomes a dependent variable (through 
the C integral). 

To study the dynamics of the reduced secular system in a 
global manner, we follow a concept of the so called representative 
plane of initial conditions introduced by [Michtchenko & Malho-[ 
[to| ( [2004| ). The representative plane (5) comprises of points in the 
phase space which lie on a specific plane crossing all phase trajec- 
tories of the system. In the cited paper, it has been shown that a 
good choice of the 5-plane flows from the condition of vanishing 
derivatives of the secular Hamiltonian over A05. It is equivalent to 
the symmetry of secular interactions with respect to fixed apsidal 
line of a selected orbit and it means that time derivatives of the con- 
jugate momenta (which can be expressed by eccentricities) must 
vanish. Indeed, in accord with the law of conservation of the total 
angular momentum, the eccentricities must reach at the same time 
maximal and minimal values along the phase trajectories, hence 
ei = 0, / = 1,2. Regarding the classic problem with Newtonian 
point-to-point interactions, this condition is satisfied when A05 = 
or A05 — 71 ([Michtchenko & Malhotra[2004 l [see also ( [Migaszewski[ 
[& Gozdziewski[[2008a| )]. If we consider the generalized problem 
of two-planet system then the concept of the representative plane 
is still valid. The small perturbations which we introduce do not 
change the dimension of the phase space and the motion remains 
on perturbed Keplerian orbits, only the form of the secular Hamil- 
tonian is modified. Moreover, we may expect that qualitative prop- 
erties of the secular system may be changed because we introduce a 
few new parameters describing physical setup of the studied system 
and interactions governing its dynamical evolution. To encompass 
both the specific values of angle A05 = 0, tt, the representative plane 
can be defined by the following set of points: 

S = {e\ cosA05 X e2\ei,e2 G [0, 1)}, 



http : //exoplanet . eu 



where A05 = %, ei cosA05 < for the left-hand half-plane (called 
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Figure 1. The long-term, secular evolution of eccentricities of the innermost planets of known non-resonant extrasolar systems. Red curves are for the 
point-mass planets, interacting mutually through Newtonian forces. The blue curves are for generalized model of motion, including additional effects (general 
relativity and the quadrupole moment of the star). For all tested planetary systems, we fix Trot = 30 days, ki = 0.02. The names of parent stars are labeling 
each panel together with a number of planets written in brackets. 



the ^Ti-plane, from hereafter), and AG5 = 0, ei cos AG5 > are for the 
right-hand half-plane (the 5o-plane). Hence, the sign of ^icosAG5 
on the X-axis tells us on the value of angle A05. 



3.2 A general view of the representative plane 

To illustrate the influence of additional perturbations on the secular 
dynamics of the classic model, we compute a set of the 5-planes 
for wide ranges of orbital and physical parameters of the studied 
two-planet configurations. Before we discuss these results, at first 
we explain how the 5-plane illustrates the dynamical structure of 
the phase space. Conveniently, Figure [2] which is derived for spe- 
cific orbital parameters (given in the caption), reveals all relevant 
features in a condensed form and its description can be regarded as 
a guide useful for the further analysis. 

We start from the left-hand panel of Fig. [2] Smooth half-ellipse 
like curves marked with thin lines are for the levels of the secular 
energy of the generalized problem. The straight lines are for the 
collision line of orbits defined through ai(l±^i)=^2(l — ^2) - The 
relative magnitude of the corrections to the secular Hamiltonian is 



represented by contour levels of the following coefficient: 



k(^1 cosAG5,^2) 



a(i^R)/aGi+a(i^M)/aGi 



(40) 



i.e., the ratio of the apsidal frequency of the inner pericenter in- 
duced by {^r-\-^m) relative to the "natural" apsidal frequency 
in the point mass Newtonian model. Regions of the 5 -plane where 
K > 0.1, are color-coded, according to the levels of constant K. Yel- 
low color (light gray in the printed paper) encodes K ^ 1, meaning 
that the apsidal frequency of the innermost pericenter forced by 
perturbing GR-i-QM corrections is larger than the relative pericen- 
ter frequency AG5 caused by the secular NG interactions between 
planets. 

The thick curves defined through: 



dGl 



0, 



(41) 



can be attributed to the stationary solutions of the reduced, two- 
dimensional system with (Gi,AG5) -variables, because coordinates 
of all points of these curves must also satisfy 



aAG3 



0, 



(42) 
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Figure 2. The representative plane {e\ cosA03,g2) for the two-planet system. The elements are mi = 1.359 mj, m2 = 0.453 mj, ai =0.1 au, ^2 = 1-0 au. 
Stellar mass is 1 Mq, rotational period is TJ-ot = 30 days, stellar radius is 1 Rq, ki = 0.02. The thick curves are for the stationary solutions: red curves 
are for stable equilibria in the classic model, blue and violet curves are for stable and unstable equilibria of the generalized model, respectively. The left 
panel: Color contours are for the ratio of the apsidal frequency induced by the general relativity and quadrupole moment to the apsidal frequency caused by 
mutual interactions between planets, K. The scale of color code is limited by K = 1, meaning that the corrections become more important than the point mass 
Newtonian gravity. In general, the yellow colour may refer to K > 1. Energy levels of the generalized model are marked with thin lines. The right panel: A 
few particular energy levels of the generalized model are labeled with a — f, accordingly. See also phase diagrams for these energy levels which are shown in 
Fig.[3]and Fig.|4] Curves of equilibria are labeled: mode I are for aligned apsides (A03 = 0), mode II are for anti-aligned apsides (A03 = n), UE means unstable 
equilibria. The thin half-ellipse like curves are for energy levels computed for the classic model (red) and generalized model (green). More details can be found 
in the text. 



according to the definition of the 5-plane. Such solutions are peri- 
odic orbits of the full secular system (Eqs.[38}j39|. We examine the 
Lyapunov stability of these equilibria with the help of the Lyapunov 
theorem, adopting i^ec as the Lyapunov function in the cases when 
they correspond to its maximum (or a minimum) in the reduced 
two-dimensional phase space. Because the investigated dynamical 
system has one degree of freedom, the stable or unstable equilibria 
can be easily identified as extrema or saddles of the secular Hamil- 
tonian in the representative plane, respectively. 

The positions of equilibria [or stationary modes fMicht chenko] 
|& Malh otra 2004)] help us to distinguish between different types 
(families) of orbits characterized by librations of AG5 around a 
particular value (libration center). The red, thick curves drawn in 
both panels of Fig. |2] are for stationary modes in the classic model 
( [Michtchenko & Malhotra|2004j ), while the equilibria in the gener- 
alized model are drawn with thick, blue curves. In the So half -plane, 
these stationary solutions are classified by [Michtchenko & Malho-| 
|tta| ([20041 as mode I equilibria, and they are characterized by libra- 
tions of AG5 around in the neighboring trajectories. In the 5% half- 
plane, we can find also Lyapunov stable mode II solutions related 
to librations of AG5 around n in close trajectories. Some parts of the 
equilibria curves are marked with violet color (for the generalized 
model). These points denote unstable equilibria (UE) accompanied 
by the true secular resonance (the TSR from hereafter). Such un- 
stable equilibria in the 5o -plane are discovered by [Michtchenko &| 
(Malhotra, (2004 ) in the secular classic coplanar model of 2-planets. 
Obviously, mode I and mode II solutions known as generic features 
of the classic model, exist also in the generalized problem. How- 
ever, their positions in the phase space may be heavily affected by 
apparently subtle GR and QM perturbations. 

Note that along the red thick curves representing equilibria in 
the classic model, K is undefined. 

Now, let us examine the right-hand panel of Fig. |2] In this 
plot, besides levels of the secular energy of the classic model (red, 
thin curves), we also plot such levels for the generalized problem 



(green, thin curves). We can observe a significant discrepancy be- 
tween the shapes of contour levels of both Hamiltonians. 

In this plot we mark also a few specific levels of O-Q^c- 
'Ea = -2.247, ^Eb = -2.24955, = -2.25, = -2.2525, 'E, = 
—2.256, "Ei = —2.26, respectively; here, we skipped all constant 
terms (the Keplerian part, indirect term, constant relativistic term) 
in the full secular Hamiltonian, and the energy values are given 
in terms of 10~^ Mq au^yr"^ when 1 Mq, 1 au, and 1 sideral year 
are taken as units of mass, distance, and time, respectively (then the 
Gauss constant k = In). Because the motion of the secular system 
must be confined to fixed energy level, points at which a particu- 
lar energy level crosses the curve representing stationary solutions, 
tell us on distinct equilibria and their bifurcations. For instance, fol- 
lowing energy level labeled with a in the right-hand panel of Fig. 2, 
we see that it crosses the equilibria curves at two points. The first 
one, found in the So half-plane, corresponds to mode I solution. 
The second cross-point is found in the S% half-plane, and it marks 
the mode II equilibrium. 

To estimate the precision of the analytical theory, the energy 
levels in Fig. [2] are calculated wit h exact semi- analytical averaging 
( [Michtchenko & Malhotra|2004| ) that makes use of precise adap- 
tive Gauss -Legendre quadratures ([Migaszewski & Gozdzie wskij 
|^08b ). For a comparison, the stationary modes are computed with 
both methods: we recall that thick curves represent equilibria cal- 
culated with the semi-numerical averaging while the thin curves 
(over-plotted on them) are for stationary solutions calculated with 
the help of the analytical theory outlined in Sect. 2 (the right-hand 
panel of Fig.[2|. The results of these methods are in excellent ac- 
cord. The most significant discrepancies between the results appear 
in the So half -plane (for AG5 = 0), in the regime of large eccentricity 
e2- In such a case, the secular series representing diverge over 
the anti-collision line (marked with red, straight line). This prob- 
lem is discussed in detail in ( [Migaszewski & Gozdziewski|2008a] ). 
Obviously, over this anti-collision line, r2 > ri at some parts of or- 
bits, breaking the underlying assumption that we require to expand 
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the inverse of the mutual distance in convergent series. Moreover, 
we demonstrate that the analytical theory reproduces the dynamics 
of hierarchical configurations up to very large 62. Another, direct 
test of the precision of the analytic theory is given in Sect. 3.4. 

3.3 Phase diagrams and the non-linear secular resonance 

To investigate more closely the secular dynamics related to the new 
mode II and UE solutions, and to understand the structure of the 
phase space in more details, we compute a number of phase plots 
(or phase diagrams). The secular energy is kept constant and we 
draw the phase trajectories in the (ei cos AG5, ei sin AG5)-, as well as 
(^2COsAG5,^2 sinAG5)-planes, choosing the initial conditions along 
the fixed energy level. The phase plots are computed with the help 
of the analytic theory. For a reference, we recall the right-hand 
panel of Fig.|2]which illustrates the 5-plane of {ei cosAG5,^2)- We 
recall that the UE of the generalized model are marked with violet 
curves. 

The phase diagrams are illustrated in Fig. [3] The top row is for 
the energy level of 'Ea = —2.247 x 10~^ which is labeled with a in 
the right-hand panel of Fig. [2] This energy level crosses the equi- 
libria curves in two points which can be recognized in the phase 
diagrams as libration centers of AG5 = (labeled with mode I), and 
of AG5 = 71 (labeled with mode II), respectively. Overall, the phase 
diagrams look like qualitatively the same as in the classic model. 
Both libration modes are separated from the circulation zone of 
AG3 by false separatrices ( Pauwels 1983, Mi chtchenko & Malhotra| 
|2004 ), i.e., the transition between each mode and the circulation of 
AG3 does not involve solutions with infinite period. This is illus- 
trated in smaller, bottom plots accompanying each phase diagram, 
which show the secular, fundamental frequency g of the mean sys- 
tem for initial conditions lying on the x axis of the respective phase 
diagram. This frequency has been determined through solutions of 
the secular equations of motion. Clearly, when the true separatrix 
crosses the ^i^2CosAG5-axis, g decreases to 0. Crossings of false 
separatrices do not lead to any discontinuity of smooth plots of g. 
The phase diagrams are much more complicated for energy levels 
labeled with b-f, which cross the equilibria curves at more than two 
points. We analyze in detail the energy level (c) which intersects the 
curve of stationary solutions in/owr points. These points can be rec- 
ognized in the phase diagrams shown in the bottom panels of Fig. [3] 
We can identify them easily in the x-axis of these diagrams (be- 
cause, such points of the 5-plane have y = ei^2 sin AG5 = 0). Starting 
at the 5o -plane and following the energy level counter-clockwise, 
we have a stable mode I equilibrium surrounded by large zone of li- 
brations of AG5 around which corresponds to single crossing point 
of the energy level and the equilibria curve in the 5o -plane. In the 
Jj^-plane, we have three such points, two of them are Lyapunov 
stable, and one point in the middle is an unstable equilibrium. This 
part of the phase space, as seen in Fig. |3] encompass figure-eight 
shaded area, involving two islands of librations (TSRs, or elliptic 
points) around AG5 = 7i, and the hyperbolic UE point lying in the 
middle between them. Both libration centers are characterized by 
AG5 = 71. Because they are related to stable equilibria separated by 
hyperbolic structure of the UE, two parts of the phase curve sur- 
rounding the islands of the TSRs and that meet in the UE, must 
form a real separatrix. The whole structure may be still surrounded 
by a zone of librations of AG5 around 7i. It is shaded in light- gray. Let 
us note, that this mode II libration area is confined to the -plane 
(hence, in this particular case, angle AG5 does not pass through 0). 

A sequence of (^2COsAG5,^2sinAG5)-diagrams for the secular 
energy levels a-f are shown in Fig.|4] Looking at these levels plot- 



ted in the 5-plane, we can now follow a development of dynamical 
structures related to the different modes of motion. In particular, 
phase diagrams Fig.|4|D and Fig.|4fi reveal bifurcations of mode II 
which emerge the UE and TSR solutions. Clearly, the bifurcations 
may be identified with points at which the given energy level is 
tangent to the equilibria curves in the 5 -plane. 

The phase diagrams assure us that the secular dynamics of the 
generalized model can be much more complex and rich due to the 
GR and QM corrections to the NG Hamiltonian than the secular 
dynamics in the classic model. 

3.4 Numerical test of the analytic secular theory 

Finally, we illustrate limitations of the secular theory and we com- 
pare its results with the outcome of the direct numerical integra- 
tions. This comparison is also directly related to the dynamical sta- 
bility of the planetary system. 

First, we constructed the numerical model of the generalized 
system independently on the analytical model. We wrote the equa- 
tions of motion with respect to the Jacobi reference frame, using 
formulation of |Mardling & Lin| ( |2002| ) who call it the direct code. 
In this code, the GR acceleration is modeled with the PPN for- 
mulae given in ( [Kidder |p^95j . In that way we have a possibility 
to check the analytic theory in completely independent way, which 
also prevents copying logical errors which could be done during the 
averaging. After some experiments, we also found that the choice 
of the reference frame (e.g., related to Jacobi, Poincare, or classic- 
astrocentric coordinates) is in fact irrelevant for the results of this 
test. We also do not account for the difference between the osculat- 
ing and the mean elements. 

Using the direct code, we integrated numerically a few phase 
diagrams for nominally the same initial conditions and parameters 
used to draw plots in Fig.|4^-f with the help of the analytic, secu- 
lar model. We computed osculating elements related to the Jacobi 
reference frame over a few secular cycles. The numerically derived 
phase curves are drawn with black, filled circles in respective pan- 
els of Fig. [4] The solutions obtained with the analytic theory are 
over-plotted on these numerical solutions with thiner, green curves. 
Subsequent five panels of Fig. |4^-e reveal that the agreement of 
both sets of solutions is excellent. The analytic theory reproduces 
qualitative features seen in the phase plots, and their structure with 
great accuracy. We find that both solutions coincide even in the 
regime of large eccentricities. Remarkably, the direct code integra- 
tions last over CPU time which is by a few orders of magnitude 
longer than the calculations carried out with the help of the analyt- 
ical theory. 

However, in the last panel of Fig. |4f we can observe signif- 
icant deviations of the analytic solutions from the numerical the- 
ory, particularly in the outer parts of the phase diagram. In fact, 
in this case ^2 is so large that the assumptions of the secular the- 
ory are broken. After examining the 5-plane (Fig. |2]), we can see 
that the energy level corresponding to the last panel passes close 
to the collision line. To illustrate the real border of the dynam- 
ical stability, we examined the dynamical character of solutions 
in the 5 -plane with the help of the Spectral Number technique 
( [Michtchenko & Ferraz-Mello|2001| ). This simple FFT-based algo- 
rithm makes it possible to distinguish between chaotic and regular 
solutions. The dynamical maps shown in Fig. [5] are constructed by 
counting the number of frequencies in the FFT- spectrum of the time 
series, {o{t) = a{t)Qxpi'ki{t)}, where ai{t) and 'ki{t) are temporal 
canonical semi-major axes and mean longitude of each planet. The 
number of peaks (the Spectral Number, SN from hereafter) in the 
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Figure 3. Phase diagrams computed for two-planet secular system described with the same parameters as used for the construction of Fig. [2] Each phase 
diagram is accompanied by a smaller plot of the fundamental frequency g of the secular solutions calculated for initial conditions lying on the _y = ^1^2 cos A03- 
axis. Panels in the top row are for the secular energy Es^ = —2.247 x 10~^, panels in the bottom row are for = —2.25 x 10~^. These levels are labeled 
in Fig.|2]with (a) and (c), respectively. The left-hand panels are for the {e\ cos ACH, ei sin Acn) -plane, the right-hand panels are for (^2 cos ACH, 62 sin ACH) -plane. 
Shaded regions mark different libration zones around mode I (A03 = 0) and mode II (A03 = tt), respectively. UE marks unstable equilibria accompanied by 
libration zones around AG3 = 71 and centered at the TSR solutions. The true separatrices are indicated by ^ ^ 0. See the text for more details. 



spectrum over some noise level tells us on the character of orbit. 
Orbits with large SN (grater than 1000) are very chaotic, while the 
SN ~ 1 means small number of frequencies and a regular, quasi- 
periodic phase trajectory. Each point in the dynamical maps repre- 
sents a phase trajectory that was integrated over ~ 10^ ^2- Although 
such time span is relevant for the short-term dynamics only, calcu- 
lations took a very long CPU time (a few days on 24 AMD-CPU 
cores). The results are shown in two panels of Fig.|5] The left panel 
is for the classic model (only NG interactions are included), while 
the right panel is constructed for the generalized model. 

Let us analyze the left-hand panel of Fig. [5] for the classic 
model. Solutions, which appear strongly chaotic are marked with 
colors (darker point means larger SN and more chaotic system). 
Clearly, the border of stable motions is irregular and is shifted to- 



wards small 62 by 0.1-0.2 with respect to the formal, geometrical 
collision line bordering the triangular region (it is drawn in both 
panels of Fig. [5}. The thick red curves mark the equilibria of mode 
I and mode II, respectively. Shaded regions are for the initial con- 
ditions in the 5-plane corresponding to orbital configurations with 
librating AG5. In the right-hand panel of Fig. [5] we show the equi- 
libria curves of the generalized model and libration zones of AG5 
associated with these equilibria. We mark again the SN signature 
of the short-term dynamics. We note that the border of stability is 
quite different from that ones of the classic model (compare with 
the left-hand panel of Fig. [5]). This is a very clear example showing 
that apparently subtle GR-i-QM effects may affect the short-term 
stability of the system in a significant way. 

The dynamical map for the generalized model (Fig.[5] the right 
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Figure 4. Phase diagrams in the {e2 cos A03, 62 sin A03)-plane drawn to illustrate a comparison of the secular evolution of the mean system (green, thin curves) 
with the numerical solutions of full, unaveraged system (larger, filled black circles). The orbital parameters are the same as in Fig.[2| Subsequent panels labeled 
with a-f correspond to relevant energy levels marked and labeled in Fig.[2| accordingly. 



hand panel) helps us to identify the source of unstable behavior 
seen in Fig.[4f, revealing that some initial conditions lead to erratic 
and irregular behavior. In the dynamical map, we mark two levels 
of i^ec corresponding to phase diagrams drawn in Fig.|4^,f, respec- 
tively. The energy level (e) lies entirely in the regular region, in 
spite that 62 may reach values as large as 0.8. The neighboring level 
(f) may touch the unstable zone, and that is why the orbital evo- 
lution at this energy level may become very unstable. Because the 
motion must be confined to the fixed energy level, due to the secular 
evolution, some initial conditions may be transported to the chaotic 
zone (by excitation of the eccentricity) during roughly a half of the 
secular period. Then the short-term, strong chaos can destabilize 
the system immediately. Following the fixed energy curve, we can 
also identify three islands of stable motions. The first one lies in the 
right half-plane of the phase diagram and is associated with libra- 
tions of AG5 around 0, see the neighborhood of the corresponding 
cross point in the 5o-plane, Fig. |2] In the Jj^-plane, we can find 
corresponding unstable equilibrium and two stable solutions with 
associated libration islands shown in Fig.|4f (see the left-hand half- 
plane of the phase diagram) around (0,-0.1) and (0,-0.78), re- 
spectively. Finally, to illustrate the development of the secular 
instability, we solved the equations of motion of the full system, 
starting very close to the UE lying on the energy level between 
levels e and f, as marked in Fig. [2] and Fig. |4^,f. The solution is 
illustrated in the phase diagram in the left-hand panel of Fig. [6] 
For a reference, the numerical solution is over-plotted on the an- 
alytically derived separatrices of the UE. The corresponding time 
evolution of the orbital elements are illustrated in the right-hand 
panels of Fig. |6] Clearly, during quite a long time the full system 



stays close to the UE, but after ~ 10^ yrs it follows a trajectory 
close to the inner separatrix, and finally begins to move close to the 
outer separatrix, approaching large 62- During this evolution, we 
observe not only very irregular behavior of a2 and both eccentrici- 
ties, but also AG5 changing from large amplitude librations around 
to circulations. Although the configuration seems bounded during 
many secular periods, such behavior may be classified as strongly 
chaotic. 



4 PARAMETRIC SURVEY OF TWO-PLANET SYSTEMS 

The characterization of the phase space with the help of the rep- 
resentative plane can be very useful to conduct a survey of the 
basic features of the secular dynamics. In particular, we want to 
understand how it depends on the physical and orbital parameters 
governing the magnitude of the GR and QM interactions. In subse- 
quent diagrams of the 5-plane, we will always mark the collision 
and anti-collision lines. In this way, we can determine the border 
of validity of the analytic approach. Yet to derive the stationary 
modes possibly exactly, in the whole permitted range of eccentric- 
ity, we compute their locations with the help of the semi- analytical 
averaging algorithm. 

4.1 Dependence of the secular dynamics on the masses 

The results of the survey of the secular dynamics of two-planet 
systems, including GR and QM interactions, for varied planetary 
masses, are illustrated in Fig. [7] We fix the system parameters as 
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Figure 5. The 5-planes for the classic model (the left-hand panel) and for the generalized model (the right-hand panel) for the two-planet system analyzed in 
Fig. [2] The thick, red curves mark stationary modes in the classic model, the thick blue curves are for the equilibria in the generalized model. Shaded areas 
indicate zones of Acn librations. The thick red lines are for the collision line of orbits. Colors code log SN of the outer orbit, characterizing solutions derived 
numerically with the help of the direct code. Black points are for strongly chaotic solutions with the logSN ~ 3, white color is for regular solutions with 
logSN ~ 0, intermediate values are marked with the color scale above the panels, accordingly (see the text for more details). 
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Figure 6. The phase diagram in the (^2 cos ACL5,e2 sin A03) -plane (the left-hand panel) drawn to illustrate a development of the secularly unstable behavior 
of the full system. The numerical solution of the full system (thin, black curve) is over-plotted on the analytical, secular solution (gray, thicker curve). The 
right-hand panels illustrate the relative changes of semi-major axes (top panel), eccentricities (middle panel) [grey curves are for outer orbit, black curves are 
for the inner orbit], and the apsidal angle Acn (bottom panel). The initial orbital parameters are the same as in Fig. [2] The energy level of this solution lies 
between levels e-f marked in Fig. [2] accordingly, compare it also with phase diagrams Fig.|4^,f. 



follows: the mass of the parent star is mo = 1 Mq, the equatorial 
radius of the star Rq = \ Rq, and Tj-ot = 30 days, ki = 0.02 (then 
J2 ~ 10~^), the semi-major axes of the planets are ai =0.1 au, and 
^2 = 1.0 au, respectively. Hence, we consider a typical hierarchical 
configuration with the orbital periods ratio ~ 30. In this experiment, 
the planetary masses are varied, but their ratio is kept constant, 
mi/m2 = 3. For a reference, the mass of the inner planet is labeled 
in the top-left corner of each respective panel in Fig.[7] Basically, 
in this test we can also analyze the effect of unknown inclination 
of the co-planar system on the long-term dynamics and stability. 
However, as we show in a recent work regarding the 14 Herculis 
planetary system |Gozdziewski et aLpOOS ), planetary masses de- 
rived from observations do not necessarily always scale according 
to the law of the mass-factor 1 / sin /. If the minimal masses are 



large then the mutual interactions in low-inclination configurations 
can strongly modify the RV signal and even the mass hierarchy 
may be reversed in the orbital fits. 

Figure [7] reveals that curves representing stationary modes, 
which are known in the classic problem, are usually significantly 
shifted and/or distorted. Also new features of the 5-plane appear 
and it can be seen in the top-left panel of Fig. [7] In general, the 
distortions of equilibria curves are the more stronger when the 
masses are smaller. It is quite straightforward to explain this effect. 
When the planetary masses decrease, also their mutual interactions 
(scaled by mim2) are decreasing. Yet the pericenter frequency in- 
duced by (i^o) is scaled by the mass product. Simultaneously, the 
GR and the spin-induced apsidal frequencies do not depend on the 
planetary masses directly (see Eq.[3T]and Eq.[35] respectively), and 
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Figure 7. The representative energy planes {e\ cos ACiJ, 62) for A03 = (the right half-planes,5o) and A03 = tt (the left half -planes, S%). Color contours are for 
the ratio of the apsidal frequency induced by the general relativity and quadrupole moment to the apsidal frequency caused by mutual interactions between 
planets. Thick curves mark positions of stationary solutions. Red lines are for stable equilibria in the classic model. Blue and violet curves mark the positions 
of stable and unstable equilibria of the generalized model. The thin lines are for the secular energy levels of the generalized model. The thick, skew lines 
are for the collision lines defined with a\{\ ±ei) = a2{l — 62). Parameters of these systems are as follows: mo = IMq, a\ = O.lau, a2 = lau, Rq = IR©, 
Trot = 30 days, ki = 0.02. Each panel was calculated for varied planetary masses, under the condition of constant mass ratio mi /m2 = 3. The mass of the inner 
planet is written in the top left corner of each panel. 



they can be regarded as approximately constant in the given mass 
range. Therefore K increases with decreasing mi and m2. Then, 
also the assumptions of the secular theory are better fulfilled. 

Some parts of the stationary curves in the Sn half-plane (for 
AG5 = 7l) comprise of unstable equilibria (they are marked with vi- 
olet color). As we mentioned already, to the best of our knowl- 
edge, such solutions are yet unknown in the literature. Similarly to 
the non-classic equilibria discovered by [Michtchenko & Malhotra] 
( |2004| ), these solutions are accompanied by the TSR solutions and 
correspond to saddles of the secular Hamiltonian. The behavior of 
neighboring solutions tells us that they are Lyapunov unstable. This 
has been analyzed in Sect. 3. 

Actually, the sequence of panels in Fig. |7] illustrates a charac- 
teristic development of curves representing the equilibria, includ- 
ing the UE solutions. When the masses are relatively large (see the 
top-left panel of Fig. [7]), the equilibria curves are distorted and the 
unstable equilibria appear at the very edge of the -plane, in the 
range of moderate and large values of eccentricity. On contrary, in 
the classic model, the UE solutions can appear (in fact, they were 
found) only for AG5 = ( [Michtchenko & Malhotra||20 04). Seem- 
ingly, the new UE branch located in the -plane is specific only 
for this model. When the masses decrease then K grows (so the 
GR+QM effects become comparable in magnitude to the NG inter- 
actions). This leads to further distortion of mode II curves and to ex- 
panding the UE part towards moderate ^1. At some point (between 
mi ~ 1.8 mj and mi ~ 1.2 mj) both stationary curves meet in a 
bifurcation point. Here, we can explain the particular choice of pa- 
rameters used to construct Fig.|2] When the masses become smaller, 
the equilibria curves separate along 62- We note that already for 



mi ^ 1.2 mj, the 5-plane is dominated by the GR-i-QM corrections. 
We recall that in the classic case, the qualitative features of the 5- 
plane do not depend on the masses individually (Micht chenko &| 
Malhotra 2004), only on their ratio in the approximation of small 
values (see also the sequence of plots in Fig. [7]). This conclusion is 
not true anymore in the realm of the generalized model. 

4.2 Dependence of the secular dynamics on semi-major axes 

In the next experiment, we investigate the dependence of the sec- 
ular dynamics of the generalized model on individual semi-major 
axes; note that the dynamics of the classic model depend only on 
their ratio, a. The results are illustrated in Fig.jS] We proceed in the 
same manner as to draw Figs.[2]and|7] We seek for stationary solu- 
tions, and we overplot the found equilibria on color-coded contour 
levels of coefficient K. The primary parameters of the tested config- 
urations are the following: mo = 1 Mq, R = 1 Rq, Trot = 30 days, 
ki = 0.02, mi = 0.4 mj, m2 = 0.2 mj. The ratio of semi-major axes 
is kept constant, a = ai/a2 = 0.1, while the individual a 1,^2 are 
varied. For a reference, the nominal value of a 1 is labeled in the top- 
left corner in each respective panel. For decreasing ai, the deriva- 
tives of ^G^^R^^M over Gi increase, hence the magnitude of 
the respective correction to the apsidal frequency grows. Moreover, 
the GR and QM induced apsidal frequency increase faster than the 
rate of the pericenter advance forced by the NG interactions. In the 
region of 5-plane painted in yellow (light gray), the GR and QM 
perturbations dominate over the NG interactions. This region ex- 
pands quickly with decreasing semi-major axis of the inner planet, 
ai. Already for ai ~ 0.4 au (which is similar to the semi-major 
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Figure 8. The survey of the secular dynamics of two-planet system, when the semi-major axes are varied. Panels are constructed in the same way, as in 
Fig. [7] Color contours are for the ratio of the apsidal frequency induced by the general relativity and quadrupole moment to the apsidal frequency caused by 
mutual interactions between planets. The thick skew line is for collisions line defined with (1 i^i) = a2{l — 62). The masses are mo = IMq, mi = 0.4mj, 
m2 = 0.2mj, respectively, the characteristic radius of the star Rq = IR©, T^ot = 30 days, ki = 0.02. Each panel is for different semi-major axes, fulfilling 
the condition of constant a = ai/a2 =0.1. The semi-major axis of the inner planet labels each respective panel. Compared to Fig.^ an additional unstable 
equilibrium for the classic model appears and it is marked with thick green lines. 



axis of Mercury in the Solar system), the apparent corrections to 
the classic model, may contribute much more to the pericenter fre- 
quency of the innermost planet than the point mass Newtonian in- 
teractions. 

The top-left panel in Fig. |8]is for ai = 2 m and a2 = 20 au, 
respectively. For these parameters, a shift of the curve of station- 
ary solutions, when compared to the ones in the classic model, is 
already significant. In the next panel (ai = 1 au, ^2 = 10 au) the 
distortion of curves representing stationary modes is even stronger. 
Moreover, the UE mode present in the classic model in the de- 
plane, cannot be found in that half-plane anymore. Simultaneously, 
new solutions appear at the very edge of the -plane, in the range 
of large ei . For ai = 0.5 au, this branch of stationary modes is even 
more extended. Starting with this value of ai, the structure of the 
5-plane with respect to the generalized model is very different from 
that ones in the classic case. For smaller ai, the curves of station- 
ary modes are still more distorted. Clearly, these distortions cannot 
be regarded as small. This result is quite unexpected, recalling that 
the semi-major axes and the planetary masses by no means are "ex- 
treme". In spite of these "typical" parameters, the secular theories 
of the classic and generalized models lead to qualitatively differ- 
ent view of the phase space. We stress again that the notion of the 
GR and QM effects as corrections (or small perturbations) to the 
secular Hamiltonian should be understood in quite a new light. 



4.3 Dependence of the secular dynamics on the stellar spin 

In the last parametric survey, we study the dependence of the sec- 
ular dynamics in the realm of the generalized model on the stel- 



lar spin (or, effectively, on the second zonal harmonic J2). Fig- 
ure [9] illustrates the 5-plane computed for following parameters: 
mo = \ Mq, mi = 1.25 mj, m2 = 0.25 mj, ai = 0.1 au, ^2 = 1.0 au, 
Ro = lRQ,kL = 0.02. The top-left panel in Fig. [9] is for the GR 
correction only and subsequent plots are for decreasing rotation pe- 
riod (generalized model) T^ot of the star (its particular values label 
the respective plots). This sequence corresponds to increasing J2. 
The top-left panel of Fig. [9] is for a spherical, non-rotating star 
(J2 = 0, no tidal bulge), the bottom-right panel is for very fast 
rotation characteristic for a young object (then ^2 = 1 x 10~^). 
Curiously, changes of the spin encompassing that range are well 
enough to emerge new families and stationary solutions which we 
described above. They are represented, as before, by thick violet 
curves drawn in the ^ti -plane. Actually, the sequence of plots sim- 
ulates variations of flattening during the lifetime of the star. Hence, 
after skipping dissipative tidal perturbations, we may conclude that 
the structure of the phase space of the secular system (and, in gen- 
eral, also its dynamical stability) may depend not only on the ob- 
served or measured orbital parameters but also on the age and the 
spin rate of the host star. 

The structure of the 5-plane and stationary modes which are 
illustrated in Fig. [9] remind us closely Fig. |7]and Fig. |8] In fact, 
as we already mentioned, the dependence of the dynamics on 
the model parameters may be described uniformly through K, see 
Eq. |40| Indeed, the increase of the stellar spin leads to increase 
of the pericenter frequency and the nominator of Eq.|40|(note that 
other terms are constant). The same behavior of K is caused by de- 
creasing masses (see the sequence of diagrams in Fig. |7] and also 
discussion in Sect. 4.1). In that case, two terms of the nominator 
of K do not change, but the denominator decreases. Finally, if the 
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Figure 9. The representative energy planes {e\ cos A03, 62) for A03 = (the right half -planes, Sq) and A03 = tt (the left half -planes, S%). Colors are for the ratio 
of the apsidal frequency induced by the general relativity and quadrupole moment to the apsidal frequency caused by mutual interactions between planets. The 
thick lines mark the stationary solutions. The thick red curves are for stable equilibria in the classic model. The blue and violet curves mark positions of stable 
and unstable equilibria in the generalized model, respectively. Thin lines are for the energy levels of a generalized model. The thick skew line is for collisions 
line defined with (1 ±^1) = a2{l — 62). Parameters of the system are the following: mo = 1.0 Mq, m\ = 1.25 mj, m2 = 0.25 mj, a\ = O.lau, a2 = 1.0 au, 
Rq^I Rq, ki = 0.02. Subsequent panels are for different rotational periods of the star, Trot which is labeled in each respective plot. For a reference, the top 
left panel is for the GR correction only. 



semi-major axes decrease in constant ratio, K also grows because 
the GR and QM-induced correction to the apsidal frequency grows 
faster thrni the NG-forced apsidal frequency. In that sense, all point- 
mass gravity corrections governed by the described above parame- 
ter changes, modify similarly the structure of the 5 -plane. 

Considering our simplified secular model, the existence of the 
branch of stationary solutions in the regime of large ei and small 62 
may seem questionable in the real planetary configurations. In that 
regime, the tidal perturbations may become very significant for the 
secular evolution. However, the position of the bifurcation point in 
the ^Ti-plane, where the two branches of equilibria curves meet, and 
which marks the extent of the branch, is shifted towards moderate 
ei ~ 0.6-0.7 when the mass ratio mi/m2 grows. This effect can be 
observed in the geometric evolution of these branches in Figs. [8] 
[t] and [9] which form a sequence of mi /m2 = 2,3,5, respectively. 
Hence, planetary configurations in that regime may really be found 
in the Nature. 



5 THE SECULAR DYNAMICS OF THE v AND SYSTEM 

Finally, we consider the generalized coplanar problem of three 
planets. This model can be still described by the Hamiltonian writ- 
ten in the general form of Eqs. ([l}j4]), and ([7}fTT]) with N = 3. 
The averaged perturbing Hamiltonian has the form of Eqs. (T5| l, 
|T6| , |T7| l, ([30]), ( [34] ) and ( [36| . To study the properties of the sec- 
ular system, we introduce the following set of angle-action vari- 
ables related to the Poincare canonical elements ( [Migaszewski &| 



|Gozdziewski|2008a| ) : 



03 : 



1' 



:G53-G5i, G[ 
:G53-G52, G2, 
:-G53, AMD = 



(43) 



G 1 + G2 + G3 , 



where Gj = Li — Gi (see Eq. 13 1. Because 03 is the cyclic angle, the 
Angular Momentum Deficit (AMD) is conserved, and the reduced 
system has two degrees of freedom. These variables make it possi- 
ble to construct the representative energy plane in a similar manner 
as for the two-planet system. Moreover, the choice of such a plane 
is not unique. One of possible definitions can be derived through a 
direct analogy to the two-planet system case. The symmetric rep- 
resentative plane can be defined as the set of phase- space points 
fulfilling the following relation: 



^0, 



1,2,3, 



with the simultaneous conditions that a/ = 0,7i. For details, see 
( [Migaszewski & Gozdziewski|2008a| ). 

A sequence of symmetric representative planes shown in 
Fig. [To| illustrates the qualitative properties of the generalized sec- 
ular model of the \) Andromedae planetary system ( |Butler et al.| 
|1999| ). This system comprises of three planets with masses and 
semi-major axes derived through the radial velocity observations: 
mo = 1.27 Mq, mi = 0.69 mj, m2 = 1.98 mj, m3 = 3.95 mj, 
ai = 0.059 au, a2 = 0.83 au, a^ = 2.51 au. The constant of the 
AMD integral (effectively, the integral of the total angular momen- 
tum) was obtained for the following eccentricities of the nominal 
system: ei = 0.029, ^2 = 0.254, ^3 = 0.242. 
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We consider the representative plane for varying age of the 
parent star, starting with approximately 30 Myr before entering the 
ZAMS. Subsequent plots are labeled by the lifetime T relative to 
this moment taken as the zero-point of the time scale, rotation pe- 
riod Trot and stellar radius (/?o) expressed in Sun's radius. For a ref- 
erence, the classic model and with the GR correction are illustrated 
in the first two top-left panels of Fig. [To| These panels are derived 
for non-rotating, non-distorted spherical star. When the star is spin- 
ning, its second zonal harmonic may be as large J2 ~ 10~^. The 
current equatorial radius of \) Andr is approximately Rq = 1.26 Rq . 
We note that the stellar radius and the spin period of the star at 
T = — 30 Myr are taken from ( [Nagasawa & Lin||2005| ), and were 
linearly interpolated over x G [—30,0] Myr. 

The results are again quite surprising. After adding the GR 
corrections to the Hamiltonian of the classic model, the overall 
view of the phase space changes significantly. The saddle of the 
secular Hamiltonian which is present in the classic model now van- 
ishes. At its place, a new maximum of the secular Hamiltonian 
appears. Moreover, at the bottom half -plane of the representative 
plane, close to the border of the permitted region of motion, two 
new saddle points appear. 

The next diagrams illustrate changes of the structure of the 
5 -plane and a development of equilibria in a sequence simulating 
time evolution of the stellar spin. At the beginning, before the star 
enters the ZAMS, the characteristic plane reveals a sharp maximum 
and a saddle in the very edge of the region of permitted motions. 
The thin curve surrounding the maximum marks the energy level 
of the nominal system. When the rotation period increases up to 
~ 8 days, the secondary extremum (the minimum) emerges in the 
place of the saddle and it persists shortly before the ZAMS stage 
and for longer rotational periods. 

Curiously, the only feature seen in the energy diagrams, which 
survives the spin variations during the whole lifetime of the star, 
and persist in the generalized model (with the GR and QM correc- 
tions) is the stable equilibrium point related to the maximum of the 
secular Hamiltonian, which is found in the range of small ei and 
moderate 62- Curiously, the nominal system appears in the energy 
level drawn with grey (green), thick line surrounding this maxi- 
mum of i^ec- We also may notice that close to this equilibrium of 
the generalized model, a saddle of the classic model appears which 
is linearly stable. 



6 CONCLUSIONS 

In this work, we consider a generalized secular theory of copla- 
nar, A/^-planet system. Extending the model analyzed in the recent 
works devoted to the secular planetary dynamics with mutual New- 
tonian point-to-point interactions (e.g., |Michtchenko & Malhotra| 
2004, Libert & Henrard 2005, Rodriguez & Gallardo 2005, Mi- 
jgaszewski & Gozdziewski^2008aJ , we consider the influence of the 
general relativity and quadrupole moment of the parent star on the 
secular dynamics of the innermost planet and stability of the whole 
planetary system. In general, these corrections to the classic model 
still do not cover all physics governing the dynamics of such sys- 
tems. In some cases (for instance, of the short-period hot- Jupiter s), 
the tidal, dissipative torques acting between the inner planet and 
the star may be significant for the orbital evolution. However, our 
main goal is rather to extend the classic model with the perturba- 
tions that are conservative and may be well modeled in the realm of 
the Hamiltonian mechanics than to build a complete, general secu- 
lar theory. Still, this approach is useful to a wide class of systems. 



when the tidal interactions may be regarded as secondary effects, 
or are acting during much longer characteristic time- scale than the 
GR and QM perturbations. In reward, for paying the price of less 
general model, we may investigate the secular dynamics in a global 
manner. 

Our analytic model follows assumptions required by the aver- 
aging theorem. Technically, the averaging has been done with the 
help of a very simple method. This algorithm relies on appropriate 
change of integration variables. It does not incorporate any clas- 
sic Fourier expansion of the perturbing function. We obtain a very 
precise analytic model of the coplanar, A/^-planet system in terms 
of the semi-major axes ratio. It can be regarded as a generalization 
of the recent analytic secular theories of the classic model inves- 
tigated in many recent pap ers (e.g., [F ord et al. 2000 , Lee & Peale| 
.2003,.Michtchenko & Malhotra|2004^. Libert & Henrard 2005 j ). On 
the other hand, our work also covers, to some extent, the global dy- 
namics of the generalized model studied in'Mardlin g & Lin] ( |2002| ); 
Nagasawa & Lin (2005 ) with the help of the Gauss/Lagrange plan- 
etary equations of motion. We stress, however, that our investiga- 
tions are devoted to more narrow class of systems (regarding the 
conservative perturbations). 

A general conclusion which can be derived on the basis of the 
generalized theory is quite unexpected. Even in a case when the or- 
bital parameters cannot be regarded as extreme, the corrections to 
the classic Hamiltonian stemming from the general relativity and 
the quadrupole moment of the star, may affect the secular dynam- 
ics dramatically. Not only the structure of the phase space of the 
secular model changes, and new branches of stationary solutions 
appear. These solutions may bifurcate within small relative ranges 
of the parameters (for instance, when the spin of the parent star is 
changing). We show that there is no simple and general recipe to 
predict the behavior of the secular system, when the perturbations 
are "switched on". The secular dynamics of the generalized model 
becomes extremely complex and rich. For some combinations of 
the system parameters, the notion of the GR and QM effects as 
corrections to the point-mass NG interactions does not seem proper 
anymore. In some cases, these effects may be more important for 
the secular dynamics than the mutual, point mass Newtonian inter- 
actions between the planets. 

We also show that these effects may be significant for the dy- 
namical stability of planetary systems both in the short-term and in 
the secular time scales. For instance, the QM generated perturba- 
tions may directly depend on the star age and its physical parame- 
ters {ki, Rq)- In turn, these effects may strongly influence the struc- 
ture of the phase space and can imply short-term, strong chaotic 
orbital evolution during a few secular periods. 

The direct tests of the analytic theory are very encouraging. 
The results justify its great accuracy. The precision of the analytic 
calculations is very important for studying the global dynamics of 
hierarchical systems. The alternative numerical approach would re- 
quire huge CPU time because the hierarchical planetary systems 
evolve during very different time scales. Then the CPU require- 
ments of the direct numerical integrations are by orders of magni- 
tude larger than those ones needed by the analytic formulae. 
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Figure 10. The secular energy levels on the symmetric representative plane for three-planet system x> Andromede. The map coordinates are x = cosoi, 
y = e2 cos 02, where Oi , 02 are (positive values of x or _y) or tt (negative values of xoxy). Grey-colored region means forbidden motions with e^ < 0. Black 
curves are for the energy levels, green, thick levels are for the energy of the nominal v Andr system. Each panel is for a different setup of the planetary system 
model. From the left-top panel: the first panel is for the classic NG model, the next panel is for the NG+GR model. Next panels are for generalized model with 
the QM corrections parameterized by the spin rate of the parent star and its lifetime x before the ZAMS stage. 



ments on the manuscript. This work is supported by the Polish Min- 
istry of Science and Education, Grant No. 1P03D-021-29. CM. is 
also supported by Nicolaus Copernicus University Grant No. 408 A. 



REFERENCES 

Adams F. C, Laughlin G., 2006a, ApJ, 649, 992 

Adams F. C, Laughlin G., 2006b, ArXiv Astrophysics e-prints 

Agol E., Steffen J., Sari R., Clarkson W., 2005, MNRAS, 359, 567 

Arnold V. I., Kozlov V. V., Neishtadt A. I., 1993, Dynamical sys- 
tems III. Mathematical aspects of classical and celestial mechan- 
ics. Encyclopaedia of mathematical sciences. Springer Verlag 

Benitez R, Gallardo T., 2008, Celestial Mechanics and Dynamical 
Astronomy, pp 41— 1- 

Brouwer D., Clemence G. M., 1961, Methods of celestial mechan- 
ics. New York: Academic Press, 1961 

Butler R. P, et al., 1999, ApJ, 526, 916 

Butler R. P, et al., 2006, ApJ, 646, 505 

Charbonneau D., et al., 2000, ApJL, 529, L45 

Ferraz-Mello S., Michtchenko T. A., Beauge C, 2006, Regular 
motions in extra- solar planetary systems. Chaotic Worlds: from 
Order to Disorder in Gravitational N-Body Dynamical Systems, 
pp 255-+ 

Fischer D. A., et al., 2001, ApJ, 551, 1107 
Ford E. B., Kozinsky B., Rasio R A., 2000, ApJ, 535, 385 
Godier S., Rozelot J.-P, 1999, A&A, 350, 310 
Gozdziewski K., Migaszewski C, Konacki M., 2008, MNRAS, 
385, 957 

lorio L., 2006, ArXiv General Relativity and Quantum Cosmol- 
ogy e-prints, gr-qc/0609112 
Kidder L. E., 1995, Phys. Rev. D, 52, 821 
Laskar J., 2008, Icarus, 196, 1 



Laskar J., Robutel P., 1995, Celestial Mechanics and Dynamical 

Astronomy, 62, 193 
Lee M. H., Peale S. J., 2003, ApJ, 592, 1201 
Libert A.-S., Henrard J., 2005, Celestial Mechanics and Dynami- 
cal Astronomy, 93, 187 
Mardling R. A., 2007, MNRAS, 382, 1768 
Mardling R. A., Lin D. N. C, 2002, ApJ, 573, 829 
Michtchenko T. A., Ferraz-Mello S., 2001, AJ, 122, 474 
Michtchenko T. A., Ferraz-Mello S., Beauge C, 2006, Icarus, 
181,555 

Michtchenko T. A., Malhotra R., 2004, Icarus, 168, 237 
Migaszewski C, Gozdziewski K., 2008a, MNRAS, 388, 789 
Migaszewski C, Gozdziewski K., 2008b, MNRAS, submitted 
Miralda-Escude J., 2002, ApJ, 564, 1019 

Morbidelli A., 2002, Modem celestial mechanics : aspects of So- 
lar system dynamics. Taylor & Francis 

Murray C. D., Dermott S. R, 2000, Solar System Dynamics. Cam- 
bridge Univ. Press 

Nagasawa M., Lin D. N. C, 2005, ApJ, 632, 1 140 

Pauwels T, 1983, Celestial Mechanics, 30, 229 

Pijpers R P, 1998, MNRAS, 297, L76 

Poincare H., 1897, Bulletin Astronomique, Serie I, 14, 53 

Richardson D. L., Kelly T. J., 1988, Celestial Mechanics, 43, 193 

Rodriguez A., Gallardo T., 2005, ApJ, 628, 1006 

Veras D., Armitage P J., 2007, ApJ, 661, 1311 

Vogt S. S., et al., 2005, ApJ, 632, 638 

Winn J. N., et al. 2005, ApJ, 631, 1215 



© 2008 RAS, MNRAS 000,[l]{T7] 



